Skip to content

Support arbitrary xgcm topologies for boundaries and masks#22

Draft
hdrake wants to merge 43 commits into
mainfrom
topology-overhaul
Draft

Support arbitrary xgcm topologies for boundaries and masks#22
hdrake wants to merge 43 commits into
mainfrom
topology-overhaul

Conversation

@hdrake

@hdrake hdrake commented Jun 25, 2026

Copy link
Copy Markdown
Owner

Summary

Generalizes regionate from single-tile MOM6 grids to arbitrary xgcm.Grid topologies — single-tile periodic, bipolar-fold (Arctic) MOM6, and genuinely multi-tile grids (ECCOv4r4 LLC90, cubed-sphere) — by driving all grid logic from the topology-aware sectionate API (plus the bipolar north-fold boundary and multi-tile padding fixes consolidated on hdrake/xgcm@dev-v1.0.0) instead of hard-coded MOM6 specifics.

(1) Boundaries from masks — boundaries.py

  • Single-tile path preserved byte-for-byte (shared _contour_to_corner_indices; existing tests unchanged).
  • New multi-tile path: per-face boundary faces are read from a topology-aware padded mask (real-neighbour halo via the grid's own boundary/face_connections); faces internal to a tile seam are dropped; surviving directed edges are stitched into closed loops by physical coincidence. A region spanning several faces yields a single boundary loop.
  • Returns a new per-corner face index f_c_list (None for single-tile).

(2) Masks from boundaries — grid_conform.py + new geometry.py

  • Antimeridian-split MultiPolygon + per-piece regionmask rasterization OR-ed together (modeled on xcryocouple.geometry; the intent is for xcryocouple to depend on regionate, not the reverse).
  • Pole-encircling boundaries (net |Δlon| ≈ 360°) are extended to the South Pole, then rasterized through the same split+OR path.
  • Coordinate names are discovered via sectionate (get_geo_corners/get_geo_centers/coord_dict) — no hard-coded geolon/geolat.

(3) Refactor

  • f_c threaded through Region/GriddedRegion/BoundedRegion, MaskRegions, and the .gr/.grs save format.
  • Fixed the genuinely-broken .gr save/load (variable-name mismatch + unguarded children); it now round-trips f_c.
  • Generalized corner-coordinate checks (arbitrary names).

(4) ECCO LLC90 worked example — examples/5_ECCO_LLC90_multiface_regions.ipynb

  • End-to-end demonstration on the native ECCOv4r4 lat-lon-cap grid: trace the Atlantic basin into grid-conforming boundary loops spanning many tiles (across rotated seams and the Arctic cap), then verify the discrete divergence theorem with ECCO's own advective heat transport (ADVx_TH/ADVy_TH, depth-integrated).
  • The check is computed two independent ways and agrees to round-off (residual 0): cell-centred convergence via xgcm's vector-aware diff_2d_vector summed over the basin mask, vs the same transport as a boundary integral via sectionate.convergent_transport(..., positive_in=region.mask) over every traced loop — the idiomatic cross-boundary-transport calculation regionate's boundary tracing exists to feed.
  • Grid geometry and a month of native temperature fluxes are loaded from the regionate/sectionate example archive on Zenodo; no NASA Earthdata login is required (the loader no longer uses earthaccess/PO.DAAC). Flux components are promoted to float64 in the loader (the archive is float32, ~12 °C m³ s⁻¹ on these O(1e8) transports, which would otherwise dominate the residual).

Dependencies / CI

  • Add explicit shapely and xgcm >= 0.9.0. Drop the now-unused earthaccess/pydap from docs/environment.yml (ECCO data now comes from Zenodo over plain HTTP).
  • Requires the topology-driven sectionate API — the face_connections- and fold-aware neighbour finding, plus the outer-lattice multi-tile corner topology and exact padded transports — consolidated on hdrake/sectionate@topology-driven-neighbors (the successor to the earlier sectionate#47). It self-reports 0.3.3, so the floor stays sectionate >= 0.3.3 and CI installs sectionate from that ref until it is released to PyPI. This line also carries: the multi-tile is_mask_inside sign fix the ECCO example relies on — on multi-tile grids the face-frame Usign/Vsign diverge from the geometric Lsign, so without it convergent_transport(positive_in=mask) returns an inverted sign on the ECCO LLC grid (the divergence theorem closes in magnitude but with flipped sign); and the multi-tile corner-topology conservation fix at 3-tile "pinwheel" junctions where a shared corner is stored on no face (sectionate#49).
  • Requires the xgcm dev-v1.0.0 line — the bipolar north-fold boundary (padding={"Y": {"fold": ...}}, formerly xgcm#711) for single-tile tripolar-fold masks/boundaries, the face_connections halo-padding fix (xgcm#712) for correct multi-tile (lat-lon-cap / cap) padding, and the bare-DataArray other_component vector-pad fix (xgcm#749) — consolidated on hdrake/xgcm@dev-v1.0.0. The floor stays xgcm >= 0.9.0 and CI installs xgcm from that ref (--force-reinstall --no-deps, so the git build wins over any PyPI xgcm an earlier step pulled in to satisfy the floor) until it is released.
  • The connected-component labeling added in (6) uses scipy (scipy.sparse.csgraph.connected_components), which is already a hard sectionate dependency — so it adds no new dependency.
  • This PR should not be released until hdrake/xgcm@dev-v1.0.0 and hdrake/sectionate@topology-driven-neighbors are published to PyPI (the remaining packaging item — the git-ref pins in pyproject.toml/CI then become plain version floors).

Tests

46 passed, 10 skipped by default (and 56 passed with REGIONATE_REALDATA_TESTS=1). New: self-contained two-face (face_connections) fixture exercising seam-spanning (one stitched loop across faces), seam-terminating (seam edge kept), interior-cell, and MaskRegions f_c threading; plus a pole-encircling regression test; and a self-contained single-tile fold fixture (test_fold_regions.py: a fold-straddling region stitches into one loop, a genuine seam boundary face is kept, the discrete divergence theorem closes to 1e-9, and an interior region is unaffected). The five real-data end-to-end checks (the MOM6 global tripolar grid, plus three ECCOv4r4 LLC90 checks — seam stitching across rotated seams and the Arctic cap, basin boundary→velocity-face conversion, and the discrete divergence theorem), plus the MOM6 Arctic-fold check (the notebook-3 region traces as one region across the fold and closes the budget), are opt-in via REGIONATE_REALDATA_TESTS=1 with the data files present locally, so the default suite (and CI) never reaches out to the network. The multi-tile convergent_transport(positive_in=mask) sign behaviour relied on by the ECCO example is covered by a regression test on the topology-driven-neighbors sectionate branch (test_left_multitile_mask_orientation).

(5) Single-tile bipolar north fold — boundaries.py

Resolves the prior follow-up. Single-tile tripolar grids that declare the xgcm north-fold boundary (padding={"Y": {"fold": ...}}, consolidated on hdrake/xgcm@dev-v1.0.0, formerly xgcm#711) now trace fold-straddling masks into a single loop. _fold_boundaries_from_mask (dispatched via _has_fold) drops fold-seam segments interior to the region — using the grid's own fold-aware halo — and stitches the surviving arcs by physical seam coincidence (a 3-D unit-sphere key, robust to the antimeridian and poles), keeping fold crossings in the coincident-seam-corner form sectionate collapses to zero-length faces. So the boundary's velocity faces reproduce the region's flux convergence exactly; the plain single-tile and multi-tile paths are unchanged.

Demonstrated in examples/3_Arctic_heat_CM4p25.ipynb: the high-Arctic cold region — previously split into two halves whose fold-normal fluxes happened to cancel — is now one region, and its advective-heat budget closes the discrete divergence theorem to round-off (boundary integral vs interior flux convergence, float64 fluxes). The example grid loader gains an opt-in fold argument.

(6) Component-based MaskRegions (code-review finding F2)

Previously MaskRegions traced every boundary loop of the whole mask and handed each resulting region the full input mask — so a region's .mask was ambiguous (it never said which loop bounded which cells). MaskRegions now decomposes the mask into topology-aware connected components and returns one MaskRegion per component:

  • connected_components(grid, mask) (boundaries.py) labels in-mask cells on the grid's own seam-aware 4-adjacency — periodic axes, the north fold, and multi-tile face_connections all supply real across-seam neighbours, via the same topology halo (_pad_center) the boundary tracer uses. A planar labeller (scipy.ndimage.label) would instead split a seam-wrapping component and merge cells adjacent only across a wall. Labeling itself uses scipy.sparse.csgraph.connected_components (already a hard sectionate dependency, so no new dependency).
  • Each component is traced from its own single-component mask, so every returned loop belongs to that component by construction. A new MaskRegion carries .mask (exactly that component's cells — unambiguous) and .boundaries (a flat list of sectionate.GriddedSection loops; a component that wraps a seam or contains holes simply has more than one). region_dict now maps a component index → MaskRegion (previously a loop index → GriddedRegion).
  • The .gr/.grs save format gains multi-loop serialization (each loop under boundaries/loop_<k>.sec/), detected and round-tripped by open_gr.

Breaking change (version bumped 0.5.5 → 0.6.0): consumers iterate region.boundaries and integrate with positive_in=region.mask per component. The real-data mask↔boundary test was rewritten from a now-tautological "the region masks tile the input mask" check (which the labeling guarantees by construction) into the substantive per-component discrete divergence theorem — a synthetic flux field's convergence over each component's cells equals the net flux through that component's boundary loops — verified on the real MOM6 tripolar and ECCO LLC90 grids. Worked-example notebooks 2/3/5 were updated to the new API and re-executed (budgets still close to round-off, residual 0).

🤖 Generated with Claude Code

hdrake and others added 8 commits June 25, 2026 13:02
- Add regionate/geometry.py with pure-shapely antimeridian helpers:
  normalize_lon, has_antimeridian_crossing, polygon_to_lon360,
  split_at_antimeridian (ported from xcryocouple).
- Rewrite get_region_boundary_grid_indices to return a 7-tuple including
  f_c, normalizing sectionate.grid_section's 4/5-tuple return and passing
  f_c through uvcoords_from_qindices. Drop the topology= kwarg.
- Rewrite mask_from_grid_boundaries to split the region polygon at the
  antimeridian and OR per-piece regionmask rasters; discover tracer-center
  lon/lat coords via sectionate (no hard-coded geolon/geolat).
- Add regionate/tests/test_geometry.py (unit + functional).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Preserve the single-tile contourpy path byte-for-byte (shared
_contour_to_corner_indices) and add a multi-tile path: per-face boundary
faces from a topology-aware padded mask (real-neighbour halo), with
internal seam faces removed and surviving faces stitched into closed loops
by physical coincidence. Returns a new f_c_list (per-corner face index;
None for single-tile).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
- GriddedRegion/MaskRegions accept ij=(i_c, j_c, f_c) and store self.f_c;
  initiate_from_boundary unpacks the new 7-tuple from
  get_region_boundary_grid_indices.
- BoundedRegion threads f_c through grid_section/uvcoords/GriddedSection.
- Replace hard-coded geolon_c/geolat_c presence checks with a
  get_geo_corners(grid) probe (arbitrary corner-coord names).
- Fix the genuinely broken .gr save/load: write/read consistent
  lons_c/lats_c/i_c/j_c names, persist+reload f_c, set lons_uv/lats_uv and
  children only when present, and guard the children directory on load.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Add explicit shapely and xgcm>=0.9.0 (both used directly). Keep the
sectionate floor at >=0.3.3 because the topology-driven branch
(MOM6-community/sectionate#47) self-reports 0.3.3; CI installs that branch
from its PR ref until it is released to PyPI.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…geometry

- Self-contained two-face (face_connections) fixture with seam-spanning,
  seam-terminating, interior-cell, and MaskRegions f_c-threading tests.
- Optional, skipped-by-default MOM6-global and ECCO/LLC checks
  (REGIONATE_REALDATA_TESTS); ECCO asserts reciprocal-or-raises.
- Export regionate.geometry helpers from the package namespace.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
A boundary whose net longitude winding is ~360 deg encircles a pole and has
no simple lon/lat polygon. Re-introduce the South-Pole extension (sectionate's
stereographic-plane orientation convention) for that case only, then rasterize
its output through the same antimeridian-split + per-piece OR path as ordinary
regions. Normal bounded regions are unaffected (net winding ~0). Add a
regression test asserting a zonal equator line encloses the southern hemisphere
for both windings.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…e#47

- Re-execute all four example notebooks against the topology-driven stack
  (new sectionate API + regionate changes) on the real MOM6 example data;
  all run cleanly (validates the high-level API including BoundedRegion).
- CLAUDE.md: document arbitrary-topology support, the f_c face index, the
  rewritten mask/boundary algorithms, new geometry.py, and corrected .gr I/O.
- installation.rst: note the sectionate#47 requirement and dev-install step.
- Verified the Sphinx/nbsphinx docs build succeeds.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
@review-notebook-app

Copy link
Copy Markdown

Check out this pull request on  ReviewNB

See visual diffs & provide feedback on Jupyter Notebooks.


Powered by ReviewNB

hdrake and others added 15 commits June 25, 2026 13:40
examples/load_example_model_grid_gfdl.py has been dead code since the
Zenodo-data migration (commit 798528c): nothing in the repo imports it or
its functions (load_CM4p25/load_OM4p5/fix_grid_coords), and it only runs
against GFDL-internal /archive paths. Removing it; git history preserves it.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
New example demonstrating regionate on the ECCO lat-lon-cap (LLC90) grid -- a
real 13-tile face_connections topology:
- examples/load_example_ECCO_grid.py: downloads the ECCO geometry from PO.DAAC
  (earthaccess), renames XC/YC/XG/YG to geolon*/geolat*, builds the LLC90 grid
  with canonical face_connections, and converts the native MITgcm 'left'
  staggering to symmetric via sectionate.gridutils.symmetrize.
- examples/5_ECCO_LLC90_multiface_regions.ipynb: traces a region straddling the
  tile-1/tile-2 seam -> a single boundary loop whose face index spans both
  tiles; plots the mask and stitched boundary. Added to the docs toctree.
- test_real_grids.py: replace the ECCO stub with a concrete seam-stitching test
  (runs when the geometry file is present; skipped otherwise).
- docs/environment.yml: add earthaccess + pandoc (nbsphinx render dep).

Requires sectionate.symmetrize (MOM6-community/sectionate#47 branch, local commit).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
The multi-tile stitcher matched arc endpoints only by lon/lat coincidence, which
fails across rotated (e.g. X->Y) tile seams: there the two tiles' symmetric
corners are offset by one cell and never coincide, so a region spanning such a
seam fragmented into one loop per non-rotated block (e.g. the Atlantic split
into east [1,2] and west [10,11] halves).

Add a topology-aware fallback: when neither the exact corner nor a coincident
corner continues the loop, join to an edge starting at a grid-ADJACENT corner
via sectionate.build_neighbor_maps. This closes rotated-seam and Arctic-cap
crossings. Single-tile and non-rotated multi-tile output is unchanged.

Update the ECCO example to trace a real Atlantic basin: its perimeter is now a
single loop spanning all five tiles it touches [1,2,6,10,11] (was fragmented).
Add a real-data regression test for the rotated tile-2/tile-10 seam.

Needs xgcm with face_connections padding fix (xgcm#713) for the cap/rotated
seams, and sectionate.symmetrize for native ECCO 'left' staggering.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Add a CI-runnable synthetic 2-tile grid with a rotated X->Y seam, built on the
MITgcm 'left' corner and converted via sectionate.gridutils.symmetrize. The two
tiles use disjoint coordinate blocks so the seam's corner representations never
coincide in lon/lat -- a region spanning the seam can only be stitched through
the grid topology (build_neighbor_maps), giving CI coverage of the rotated-seam
join that previously was exercised only by the (data-gated) ECCO test.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…titch

Replace the cell-based multi-tile tracer with the planned decomposition, which
isolates the only bespoke step (cross-seam stitching) and produces boundaries
that are grid-adjacent everywhere -- so they convert to velocity (u,v) faces and
obey the discrete flux-divergence theorem across rotated seams AND the Arctic cap.

Algorithm (regionate/boundaries.py):
  1. contourpy traces each face's mask boundary (reusing the single-tile path).
  2. a segment is internal (cut) iff BOTH cells it separates are in-mask in the
     topology-aware padded mask -- the entire seam-cutting rule.
  3. surviving open arcs are stitched by the global cell-set at their endpoints
     (cells are globally unique, so no coordinate coincidence is needed).
  4. face-local corners are converted to native (f,j,i) by a cell-set-matched
     adjacency walk; seam crossings the cell-set over-merges are repaired by
     inserting the corner where the two tiles meet; the loop is closed once.
Works for 'left' (MITgcm/ECCO), 'outer' (symmetric) and 'right' grids via
sectionate's corner_position/corner_offset; single-tile uses the same offset fix.

ECCO example now loads the native 'left' LLC90 grid directly (no symmetrize) and
uses the published regionmask North+South Atlantic basin (geographically correct,
no Pacific contamination). The basin traces to one 720-corner boundary across
six tiles and converts to 720 velocity faces. New real-data test asserts the
Atlantic boundary converts to velocity faces (the transport-consistency check).

Requires the topology-driven sectionate (MOM6-community/sectionate#47): native
'left' support + corner_position/corner_offset; and xgcm with the face_connections
padding fix (xgcm#713) for the cap.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…atitude-circle section

- atlantic_basin_mask: union the Atlantic-sector Natural Earth basins (the open
  Atlantic is split into several, notably the Sargasso Sea filling the central
  North Atlantic -- selecting only '*Atlantic*' left a large hole). Fills the
  basin, follows coastlines, excludes the Pacific.
- Actually verify the discrete divergence theorem (do not just assert it): on a
  region crossing a non-rotated tile seam, the flux convergence summed over the
  cells equals the net flux through the traced boundary's velocity faces to
  round-off (residual ~1e-14). Added as a regression test, and demonstrated in
  the notebook (with an honest note that a cell-centred divergence across the
  LLC's rotated seams is subtle -- u rotates into v -- and is handled by
  sectionate.convergent_transport).
- Notebook: distinct styling for boundary vs interior points; add a 26N
  latitude-circle section (curve='latitude circle') contrasted with great circle.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…sion pins

- examples/5 (ECCO LLC90): trace the full Atlantic basin, draw ALL boundary loops
  (so every masked cell is enclosed), and verify the discrete divergence theorem
  over the WHOLE basin (residual ~1e-14) using xgcm's vector-aware diff_2d_vector
  with an arbitrary transport field; remove the latitude-circle section.
- boundaries.py: call the module-level xgcm.padding.pad directly, dropping the
  dead hasattr(grid,"pad") fallback (the pinned xgcm exposes only that API).
- tests: add a real-ECCO Atlantic-basin divergence-theorem regression; compute
  cell-centred convergence with grid.diff_2d_vector (correct across rotated seams,
  where a per-component scalar grid.diff would be wrong -- expected xgcm behaviour).
- deps (CI git refs + pyproject/meta comments): require the xgcm face_connections
  padding fix (hdrake/xgcm@fix-faceconnection-pad-712) and the topology-driven
  sectionate API + 'left'/rotated-seam corner-neighbour fix (sectionate#47), both
  needed for correct multi-tile (lat-lon-cap / cap) results.
- docs: re-execute all example notebooks.

The rotated-seam boundary tracing bug this exercises is fixed upstream in
sectionate#47 (build_neighbor_maps corner-neighbour correction on 'left'/'right'
grids); validated to machine precision on the Atlantic basin and the Arctic cap.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01E9R9Y9HxT5ZQyq6hKx4Hjm
Load the ECCO LLC90 geometry and a month of native 3-D temperature
fluxes from the regionate/sectionate Zenodo archive
(10.5281/zenodo.21051424) instead of streaming from PO.DAAC via
earthaccess -- no NASA Earthdata login required. Drop the now-unused
earthaccess/pydap deps from docs/environment.yml.

Rework the divergence-theorem check to use ECCO's own advective heat
transport: ADVx_TH/ADVy_TH depth-integrated with .sum("k"). Compare the
cell-centred flux convergence (xgcm vector-aware diff_2d_vector, summed
over the basin mask) against the same transport as a boundary integral
via sectionate.convergent_transport over every traced loop. The two
agree to round-off (residual 0). The flux components are promoted to
float64 in the loader (the archive is float32, ~12 degC m3 s-1 on these
O(1e8) transports, which would otherwise dominate the residual).

Requires the is_mask_inside multi-tile sign fix in sectionate#47.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
- test_real_grids: the ECCO end-to-end tests were effectively ungated --
  the `skipif` marker decorated the `_load_ecco` helper, where pytest
  ignores it, so the three ECCO tests ran unconditionally and triggered
  live Zenodo + Natural-Earth downloads in CI. Gate each test (and the
  MOM6 test's peer) on REGIONATE_REALDATA_TESTS=1 + a local data file, so
  the default suite never reaches the network. Default: 16 passed,
  4 skipped; with REGIONATE_REALDATA_TESTS=1 and data present, all pass.
- download_ECCO_data: stream to a `.part` temp file and os.replace on
  success (clean up on error), so an interrupted transfer can't leave a
  truncated file that a later run treats as complete; makedirs(data_dir).
- region.py: fix to_gr docstring (writes region.nc, not boundary.nc; the
  method is to_gr, not save).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…loop

regionate now traces masks on single-tile tripolar grids that carry an xgcm
north-fold boundary (xgcm#711, boundary={"Y": {"fold": ...}}). A region
straddling the Arctic fold is traced into a single boundary loop with the fold
treated as interior, rather than splitting into two halves whose fold-normal
fluxes cancel.

- boundaries.py: _has_fold + _fold_boundaries_from_mask. It drops fold-seam
  segments interior to the region (via the grid's fold-aware halo) and stitches
  the surviving arcs into one loop by physical seam coincidence (a 3D unit-sphere
  key, robust to the antimeridian and poles), keeping fold crossings in the
  coincident-seam-corner form sectionate collapses to zero-length faces. Dispatch
  prefers multi-tile, then fold, then plain single-tile; the plain and multi-tile
  paths are unchanged.
- tests: self-contained fold fixture (straddling -> one loop; a genuine seam
  boundary face kept; divergence theorem to 1e-9; interior region unaffected) and
  a gated real-grid Arctic-fold check.
- examples: load_example_model_grid gains a `fold` opt-in; notebook 3 traces the
  Arctic as a single region and closes the discrete divergence theorem to
  round-off (float64 fluxes; boundary integral vs interior flux convergence).
- deps/docs: CI and pyproject pin xgcm to tripolar-north-fold (#711, a superset
  of the merged #712 padding fix); CLAUDE.md documents the fold path.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Import _is_fold_boundary defensively so `import regionate` still works against
an xgcm without the north-fold branch (xgcm#711); such an xgcm cannot build a
fold grid anyway, so _has_fold is correctly False there.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
The tripolar-north-fold branch's setuptools_scm version is 0.1.devN, below the
xgcm>=0.9.0 floors in sectionate#47 and regionate; installing it first meant pip
clobbered it with a PyPI xgcm 0.9.0 (no fold) while resolving those. Install it
last, forced and without deps, so the fold xgcm is the final one.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
xgcm deprecates diff_2d_vector/interp_2d_vector in favor of grid.diff/interp with
other_component=. Notebook 5 (ECCO) and the ECCO divergence-theorem test now use
grid.diff({axis: component}, axis, other_component={...}, to='center'), which
reproduces diff_2d_vector exactly (ECCO residual stays 0.00e+00). The component must
be passed as a single-entry dict so xgcm carries the axis-component context for the
seam rotation. (Notebook 3's fold budget uses scalar grid.diff, correct on the
symmetric 'outer' grid where the seam flux is stored and no U<->V rotation occurs.)

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…ce_connections)

Collapse the three mask->boundary tracers into a shared front-end plus two stitch
back-ends, so a region wrapping ANY seam the grid declares is traced into one
seam-consistent boundary instead of split with spurious seam faces.

- `_trace_and_drop(grid, mask)`: the shared front-end -- contourpy-trace per face
  (nf=1 for single-tile) + drop segments interior to a seam via the topology-aware
  halo `_pad_center`. Reused by both back-ends.
- Single-tile back-end (`_single_tile_boundaries_from_mask`, generalized from the fold
  tracer): stitch arcs by physical seam coincidence -> f_c=None. Now handles walls,
  periodic-X, and the bipolar fold uniformly; the old zero-pad plain path and the
  `_has_fold` special-case are gone.
- Multi-tile back-end unchanged in logic (cell-set + neighbour-map stitch, f_c), now
  consuming the shared front-end.
- `_pad_center` fix: coerce every non-periodic/non-fold boundary to fill+NaN. `extend`
  replicates the edge cell, which would make a domain wall look like an in-mask seam and
  wrongly drop it (e.g. the southern wall circle of an inverse zonal strip). No-op for
  multi-tile (already fill) and fold (periodic+fold kept).

Behavior: non-seam-touching regions are byte-identical to before (interior box
unchanged). A periodic-X wrapping band now traces as one loop; a globe-spanning zonal
strip as two latitude circles (each closed across the seam, no radial cut) rather than
one loop cut at the seam. Budgets close as before (a property of the full loop set).

Tests: `test_gridded_regions.py` strip/inverse expectations updated to the seam-
consistent circles (order-insensitive); new `test_periodic_regions.py` (wrapping band ->
one loop, strip -> two circles, both closing the discrete divergence theorem). Fold,
multi-tile, and gated real-data (ECCO + MOM6 fold) checks unchanged and green.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
hdrake and others added 5 commits July 3, 2026 13:57
Confirms every worked example still runs and closes its budget after the boundary-
tracing unification: thickness (NB1) and heat-convergence (NB2) divergence checks
True, Arctic fold (NB3) residual 0, ECCO multi-tile (NB5) residual 0, bounded-by-
sections (NB4) clean. Regions in these notebooks (mid-latitude, Arctic fold, ECCO
Atlantic) are unaffected in value; this just refreshes the committed outputs.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…6 & ECCO)

A new worked example (notebook 6) illustrating regionate's grid-conforming geometry
for deliberately awkward idealized regions, side by side on the MOM6 tripolar
single-tile grid and the ECCOv4 LLC90 lat-lon-cap multi-tile grid.

mask -> boundary: a North Pole cap (one loop across the fold / Arctic cap); a southern
cap that is an annulus (two loops -- neither grid has ocean at -90); compact boxes that
do and do not cross the antimeridian (one loop each); a zonal band crossing all
longitudes (two latitude circles, seam stitched not cut); and nested holes -- a lake on
an island inside an inland sea (four nested loops).

boundary -> mask: a polygon crossing the antimeridian and a diamond (staircase
rasterization); shown on MOM6, with the ECCO panel noting that polygon rasterization
(regionmask) needs a 1D/2D grid and the multi-tile workaround.

Closing notes document the pole/domain caveat, seam stitching, multi-loop-per-region
(regionate#23), and the boundary->mask multi-tile limitation. Added to the docs toctree
(conf.py auto-copies examples/*.ipynb into docs/source/examples).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Replace the cell-set + BFS-bridge multi-tile back-end with a stitch on
sectionate's outer_topology: every traced corner (contoured in the outer-
lattice frame) resolves to a physical corner node, so rotated/reversed
seams, 4-face cube-vertex junctions, the pole, and grid cuts/folds are all
ordinary nodes. Edges traced from both sides of an undeclared cut or
boundary fold (LLC90's south rows under Antarctica) annihilate as interior;
the surviving directed segments chain into closed loops by node identity,
and each node is emitted as the single native corner that stores it.

Boundaries through corner points stored on no face now raise a clear error
instead of emitting fabricated indices; masks in the synthetic rotated-seam
test stay clear of its open walls accordingly, and a new test pins the
error. Seam corners are no longer duplicated in loops, so every boundary
edge is a real velocity face (updated counts in the seam-spanning test).

Real-data ECCO tests now compute the budget area side with
outer_topology(grid).padded_transports (exact: global convergence
telescopes to zero bit-for-bit) instead of grid.diff(other_component=),
whose vector halos pick the wrong component across rotated LLC seams, and
add exact mask<->boundary closure cases for the south-pole cap (boundary
fold + grid cut), the Arctic cap, nested annuli around both flavours of
4-face cube-vertex junction, and a lat-lon<->rotated seam box.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…disk on LLC90

Notebook 5's area-side convergence now uses sectionate's
outer_topology(grid).padded_transports (topology-exact; the divergence-theorem
residual on the real Atlantic advective heat budget is now exactly zero)
instead of grid.diff(other_component=), whose vector halos are wrong across
rotated LLC seams. Notebook 6's southern-cap case documents the corrected
behavior: MOM6's open southern edge still gives an annulus (two loops), while
LLC90 reaches -90 with a folded southern boundary that the tracer recognises
as interior, so the same cap traces as a single-loop disk containing the pole.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…tebooks

Point CI at the sectionate outer-lattice PR (MOM6-community/sectionate#48,
refs/pull/48/head) for the most complete multi-tile boundary treatment, and at
the combined xgcm branch tripolar-north-fold-vectorpad (bipolar north-fold
#711/#712 + the bare-DataArray other_component vector-pad fix #749). Re-executes
all six example notebooks against that environment: notebook 5's ECCO LLC90
divergence-theorem residual is exactly 0.00e+00.

Also nuances the padded-transports wording in CLAUDE.md and the test_real_grids
docstring: xgcm#749 fixes the bare-DataArray grid.diff(other_component=) path
across rotated/reversed seams (the dict form was always exact), but
padded_transports remains the exact, xgcm-independent path because it resolves
edges stored on no face (walls, the lon=-115 cut, the 4th Arctic vertex) that no
halo pad can supply a value for.

Tested: pytest with REGIONATE_REALDATA_TESTS=1 -> 35 passed.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
@hdrake

hdrake commented Jul 4, 2026

Copy link
Copy Markdown
Owner Author

Rebased this PR's environment onto the outer-lattice stack for the most complete multi-tile boundary treatment:

  • sectionateMOM6-community/sectionate#48 (outer corner-node topology + exact padded_transports; builds on #47), installed via refs/pull/48/head.
  • xgcmhdrake/xgcm@tripolar-north-fold-vectorpad = the bipolar north-fold branch (xgcm#711/#712) with the bare-DataArray other_component vector-pad fix (xgcm#749) cherry-picked on top.

Validation against that environment:

  • pytest with REGIONATE_REALDATA_TESTS=135 passed (includes the ECCO LLC90 discrete-divergence-theorem checks).
  • All six example notebooks re-executed cleanly; notebook 5's ECCO LLC90 divergence-theorem residual is exactly 0.00e+00.

Docstring/CLAUDE.md wording for padded_transports was nuanced: xgcm#749 fixes the bare-DataArray grid.diff(other_component=) path across rotated/reversed seams (the dict form was always exact), but padded_transports remains the exact, xgcm-independent path because it resolves edges stored on no face (walls, the lon=−115 cut, the 4th Arctic vertex) that no halo pad can supply a value for.

hdrake and others added 15 commits July 13, 2026 09:20
Adopt the AI Usage Policy drafted for the sibling xgcm project
(disclose AI assistance; be able to explain the diff) and add five
engineering norms adapted to regionate.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…bors

Point CI (and the version-floor comments in pyproject.toml / conda/meta.yaml /
CLAUDE.md) at the consolidated upstream dev branches:
  - xgcm  -> hdrake/xgcm@dev-v1.0.0            (north-fold #711 + face_connections
             pad #712 + vector-pad #749, all in one branch)
  - sectionate -> hdrake/sectionate@topology-driven-neighbors (223cb76; formerly #47/#48)

dev-v1.0.0 removes two long-deprecated xgcm.Grid kwargs, so migrate to the new API:
  - Grid(boundary=...)      -> Grid(padding=...)
  - Grid(periodic=False)    -> Grid(padding="fill")   (ECCO loader)
  - axis.boundary attribute -> axis.padding           (regionate/boundaries.py)
This is a pure rename with no semantic change (fill_value / face_connections / the
{"fold": ...} dict structure are unchanged).

Re-ran the full suite against the new deps (35 passed, incl. REGIONATE_REALDATA_TESTS=1
MOM6 + ECCO checks) and re-executed all six example notebooks (all error-free; the
notebook 3 fold budget and notebook 5 ECCO multi-tile budgets close to round-off).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
regionmask only accepts 1D/2D lon/lat, so building a mask from a lon/lat
boundary polygon (GriddedRegion(name, lons, lats, grid)) previously raised
"1D or 2D data required - found 3 dimensions" on a multi-tile lat-lon-cap grid,
whose center coordinates carry an extra (face) dimension.

Add regionate.grid_conform.rasterize_per_tile(grid, per_slice): it applies a
2D (lon,lat)->bool rasterization once per face and stitches the per-face masks
back along the face dimension (single-tile grids call it once). mask_from_grid_
boundaries now routes its per-piece regionmask rasterization through it, so
boundary->mask works on lat-lon-cap / cubed-sphere grids as well as single-tile
ones. The expected per-tile "No gridpoint belongs to any region" warnings are
suppressed while stitching.

examples/load_example_ECCO_grid.py: atlantic_basin_mask now reuses
rasterize_per_tile instead of hand-rolling the per-tile loop (output is
bitwise-identical: 11284 cells).

Notebook 6 cases 7 & 8 now render the ECCO panel (2572 / 590 cells enclosed)
instead of the "needs a 1D/2D grid" fallback. Tests: 36 passed (incl.
REGIONATE_REALDATA_TESTS=1), with a new multi-tile boundary->mask test; both
affected notebooks (5, 6) re-executed error-free.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Snapping a polygon boundary onto the grid now follows constant-latitude arcs
(marching in longitude) by default instead of great-circle geodesics, so a box
or basin edge between two same-latitude vertices stays on that latitude rather
than bowing poleward. Thread a `curve` parameter (default "latitude circle")
through GriddedRegion, initiate_from_boundary, get_region_boundary_grid_indices,
and BoundedRegion (parent + child sections); pass "great circle" to opt back in.

Relies on the sectionate zero-length-closure fix (topology-driven-neighbors
af6e878) so a globe-encircling boundary's 360->0 loop closure is not rejected.

Adds test_curve_default_is_latitude_circle (default == latitude circle, differs
from great circle, which bows poleward to enclose more cells). Full suite: 37
passed (incl. REGIONATE_REALDATA_TESTS=1). Re-executed the boundary-constructing
notebooks (1, 2, 4, 6) error-free under the new default.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
- F5 (integrate.py): check_global_coverage now raises (was a no-op Valueus
  never raised) and reports gap vs overlap counts instead of a leaked loop var.
- F6 (regions.py): GriddedRegions validates member types with a specific
  TypeError instead of a bare except that swallowed the real cause; drops a
  duplicate super().__init__.
- F8 (grid_conform.py, region.py): remove the dormant `along_boundary` flag
  that silently disabled pole handling; the ij+mask=None path is now pole-aware.
- F1/F9: add test_cube_rotated_seam_obeys_divergence_theorem on a new
  physically-valid cubed-sphere fixture (regionate/tests/cube_grid.py, adapted
  from sectionate's cube test). It checks the discrete divergence theorem across
  rotated seams against an INDEPENDENT xgcm.diff ground truth -- the coverage
  that previously ran only in the gated real-ECCO tests. (The existing
  rotated_two_tile_grid is stitching-only: its non-physical offset coords make
  no divergence test close on it.)
- Add test_regions_validation.py covering F5 + F6.

Investigation note: the apparent xgcm-vs-sectionate rotated-seam disagreement was
a fixture artifact of rotated_two_tile_grid; on the valid cube and on real ECCO
(105297/105300 cells) they agree. The 3 ECCO exceptions are the Arctic cube-vertex
edges stored on no face, which xgcm's halo cannot supply and padded_transports
handles -- a known structural limit, not a bug.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
- F4 (grid_conform.py): _pole_enclosing_polygon now picks the enclosed pole from
  the boundary's mean latitude, so a north-polar-cap boundary encloses the NORTH
  cap instead of its complement (equatorial/ambiguous still defaults to south).
  Adds test_pole_cap_boundary_encloses_the_hugged_pole.
- Add test_cube_vertex_junction_closes_for_nondivergent_flow: exercises the
  edges-stored-on-no-face points (the 2 unstored cube vertices + a stored one)
  that the xgcm.diff seam test cannot check (xgcm.diff isn't a valid oracle
  there). Uses the div-free-streamfunction oracle (closed-loop net transport == 0,
  valid even at the vertices). Regions wrapping all three junctions close to 0.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
roll_boundary_to_align_with_overlap rolled and gap-trimmed lons_c/lats_c but left
i_c/j_c/f_c untouched, so after align_boundaries_with_overlap_sections a region's
corner indices no longer pointed at its coordinates. Co-roll and co-trim i_c/j_c
(and f_c on multi-tile regions) in lockstep, handling the closed (first-corner-
repeated) storage of the index arrays vs the open coordinate arrays. Adds a
regression test on two adjacent regions sharing an edge.

Note: overlap matching stays on physical lon/lat coincidence, which is correct on
real grids (coincident lon/lat is the same physical corner).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
open_gr rebuilt each child from its stored lons_c/lats_c only, discarding the
i_c/j_c/f_c that to_gr writes -- so a reloaded region's children were bare
sec.Section objects, not the sec.GriddedSection instances BoundedRegion produces.
Reconstruct them as sec.GriddedSection(section, grid, i_c, j_c, f_c), mirroring
BoundedRegion. Adds a round-trip test.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Add connected_components(grid, mask): labels a boolean cell mask into
seam-aware 4-connected components, reusing the boundary tracer's own
topology halo (_pad_center) so periodic axes, the north fold, and
multi-tile face_connections all supply real across-seam neighbours. A
planar labeller (scipy.ndimage.label) would wrongly split a
seam-wrapping component and merge cells adjacent only across a wall.

This is the foundation for rebuilding MaskRegions so each region is one
connected component with an unambiguous mask (F2), replacing today's
"every region gets the full input mask".

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Rewrite MaskRegions so each region is one topology-aware connected
component, not one boundary loop. Previously every region was handed the
*full* input mask, so it was ambiguous which loop bounded which cells.

Now MaskRegions labels the mask into seam-aware components
(connected_components), and for each component re-traces its own
single-component mask so every returned loop belongs to that component by
construction. Each component becomes a new MaskRegion carrying:
  - .mask     : exactly that component's cells (unambiguous)
  - .boundaries: the list of sectionate.GriddedSection loops enclosing it
                 (a seam-wrapping or holey component simply has >1 loop).

Use scipy's connected_components for the labeling (scipy is already a hard
dependency of sectionate, so no new dependency); ~17x faster than a
Python union-find on a whole-domain mask.

Consumers that treated each loop as a separate region are updated to
iterate a component's .boundaries (test_gridded_regions,
test_multitile_regions). GriddedRegions accepts MaskRegion members.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
MaskRegion.to_gr writes the component's mask to region.nc (tagged
kind="MaskRegion") plus one gridded-section per boundary loop under
boundaries/loop_<k>.sec/, since a component can have several loops (a
GriddedRegion has exactly one). open_gr detects the kind tag / boundaries
dir and rebuilds a MaskRegion, its loops reloaded as sectionate
GriddedSections carrying their stored i_c/j_c/f_c. Round-trip test on a
seam-wrapping strip (one component, two loops).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…nsumers

Rewrite test_mom6_global_box_mask_boundary_consistency: the old assertion
(union of region masks == input mask) is now tautological, since the
labeling partitions the input mask by construction. Replace it with the
substantive check -- a synthetic flux field's convergence over each
component's cells equals the net flux through that component's own
boundary loops (the discrete divergence theorem), which actually
exercises that the traced boundaries enclose those cells.

Update the remaining realdata consumers to the component/boundaries API:
the MOM6 fold budget test and the ECCO seam/basin/divergence tests now
iterate r.boundaries (each a GriddedSection) and integrate with
positive_in=r.mask (the component's own cells) rather than the shared
full mask. All 10 realdata tests pass on real MOM6 tripolar and ECCO
LLC90 grids.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
… API; bump 0.6.0

Notebooks 2/3/5 consumed MaskRegions via the old per-loop API
(region.i_c / region.lons_c / positive_in=regions.mask, one region per
loop). Update them to the new model -- iterate each connected region's
.boundaries (each a GriddedSection) and integrate with
positive_in=region.mask -- and re-execute headless on the regionate-pr22
kernel:

  - nb2 (MOM6 advective heat): plot/integrate over every component's
    boundary loops; region-count print now reports connected regions and
    total loops. Both divergence-theorem checks pass.
  - nb3 (MOM6 Arctic fold): the fold-spanning cold region is one
    connected component; keep the mask-background plots on the full
    ds.mask (which the old code got implicitly, since every region used
    to be handed the full mask), so the boundary-vs-interior heat budget
    closes to 0.00 TW residual.
  - nb5 (ECCO LLC90): keep the component dict, flatten every region's
    .boundaries into the ranked loop list the figure draws; the basin
    divergence-theorem residual is 0.00.

Notebook 6 uses only grid_boundaries_from_mask + GriddedRegion and is
unaffected. Bump version 0.5.5 -> 0.6.0 for the breaking MaskRegions API.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Example 8 (boundary->mask of a rotated square) rasterized to two
diagonally-offset blocks rather than a diamond: GriddedRegion defaults to
curve="latitude circle", which snaps each polygon edge along constant
latitude — correct for axis-aligned boxes but degenerate for the
diamond's diagonal edges (each edge collapses onto a latitude line,
filling the bounding box's corners). Give plot_b2m a `curve` argument and
trace the diamond with curve="great circle" so its diagonal edges follow
geodesics; the enclosed mask is now the expected filled diamond and its
re-traced boundary the staircase of cell corners (verified on both the
MOM6 tripolar and ECCO LLC90 grids: 9196 and 590 cells).

Also refresh a note that F2 made stale: MaskRegions no longer makes one
GriddedRegion per loop — each connected component is one MaskRegion whose
.boundaries hold its loops.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Reverses the earlier default (9d25b99). Snapping a polygon boundary onto
the grid now follows great-circle geodesics by default, so arbitrary --
including diagonal -- edges are traced faithfully. The latitude-circle
default was a footgun: it marches each edge along constant latitude, which
is fine for an axis-aligned box but collapses a diagonal edge onto a
latitude line -- e.g. the notebook-6 "diamond" corner case rasterized to
two offset blocks instead of a diamond.

Flip the `curve` default to "great circle" in GriddedRegion,
initiate_from_boundary, get_region_boundary_grid_indices, and
BoundedRegion; "latitude circle" remains available for axis-aligned boxes.
Update docstrings and flip test_curve_default_is_* accordingly. Full suite
green (56 passed incl. REGIONATE_REALDATA_TESTS=1).

Drop the per-example curve override added for the diamond and re-execute
the boundary-constructing notebooks (1, 2, 4, 6) under the new default;
their divergence-theorem / budget checks still close (the diamond now
rasterizes correctly with no override).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant